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preface 

This  presentation  is  an  application  of  optimal  control 
theory  >to  determine  approximations  of  the  angle~of -attack 
history  that  maximizes  the  range  of  a  typical  air-to-surface 
missile  trajectory.  I  attempted  to  approach,  the  problem  as 
simply  and  as  practically  as'  possible,  thus  hoping  to  en¬ 
hance.  further  the  use  of  optimal  control  techniques  by  the 

practical  engineering  world. 

I  wish  to  express  my  appreciation  to  my  sponsor,  James  ; 

o 

P.  McCarthy,  Aerospace.  Engineer,  Aeronautical  Systems 
Division  and  my  advisor,  It,  Col.  R.  A.  Hannen,  for  their 
guidance  and  helpful  suggestions  in  preparing  this  paper. 

I ‘would' 'also  like  to  express  my  appreciation  to  Capt. 

Thomas  Ev  Mori  arty;,  Guidance  and  Control  Engineer,  Aeronau¬ 
tical  Systems  Division,  for  his  efforts  in  getting  me  inter¬ 
ested.  in  the  problem,  and  my  ex-neighbor,  Lt.  Steve  Faught, 
formerly  of  the  Digital  -Computation  Division  of  the  Aeronau¬ 
tical  Systems  Division,  for  his  help  in  writing  the  "Cal- 
compM  computer  program  used  to  make  the  graphs  in  this 

presentation. 
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Abstract 

Pontryagin's  Maximum  Principle,  coupled  with  the  con¬ 
jugate  gradient  iterative  technique,  is  employed  in  deter¬ 
mining  estimates  of  the  two-dimensional,  maximum  range 

trajectory  of  an  air-to- surface  missile.  Angle  of  attack 

is  used  as  the  control  parameter. 

The  motion  of  the  vehicle  is  described  by  four  state 
equations  including  standard  atmospheric  data,  and  lift  and. 
drag  data  obtained  from  wind  tunnel  test.  In  the  adjoint 
equations  Lagrangian  differentiation  formulas  are  used  to 
approximate  <the  derivatives  of  lift  and  drag  with  respect 
to  velocity  and  altitude. 

Two  quadratic  cost  functions  are  investigated — one 
involving-  a  linear  range  term  and  the'  other  a  quadratic 
term.  Both  include  a  quadratic  penalty  function  involving 

a  weighting  function  and  thfe  square  of  the  -control. 
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RANGE  MAXIMIZATION  OF  AN  AIR-TO- SURFACE  MISSILE 
I  .•  Introduction 

Maximizing  the  range  cf  currently  operating  ^air-.td'-.- 

surface  missiles,  thus  giving  the  launch  vehicle  more 
escape  time  and/or  distance,  is  one  problem  -that  is  of 
particular  interest  to  the  U.S.A.F.  Despite  the  current 
activity  in  applying  optimal  control  techniques  tc  missile 
trajectory  problems,,  little  ,has  been  done  in  applying  such 
techniques  to  practical  air-to-surf ace  missile  problems. 

This  paper  treats  the  use  of  such  a  technique  in  maximizing 
the  range  of  a  typical,  two  dimensional,  aiivto-surface 
missile  trajectory. 

The  conjugate  gradient  method,  is  the  iterative  tech¬ 
nique  used  in  this  investigation.  It  is  chosen  because  of 

the  speed  in-  which  it  converges,  and  the  relative  ease  in 
setting  up  the  problem  and  including  control  constraints. 

The  major  shortcoming  of  this  method  is  that  either  the. 
final  time  -must  be  known  or  time  must  be  treated  as  a  state 
variable  and  another  monotonely  increasing  variable,,  whose 
final  value  is.  known,  used  as  the  independent  variable  in 
the  state  and  adjoint  state  equations.  In  this  problem 
time  and  range  are  the  only  monotonely  increasing  variables, 
and  neither  end  condition  is  known.  Therefore,  time  is  used 

as  the  independent  variable  .and  "educated  guesses"  are  made 


GGC/EE/?0t5 


on  the  final  time.  These  '‘educated,  guesses"  are  actually 
obtained  by  first  simulating  trajectories  of  constant  con¬ 
trol  histories  and  then  making  corrections  as  more  informa¬ 
tion.  is  obtained  from  analysis. 

The  control  variable  is  the  angle  of  attack  which  de¬ 
termines  the  thrust  vector  and  the  aerodynamic  forces  act¬ 
ing  on  the  missile,  thus  determining  the  range  of  the  tra¬ 
jectory.  The  angle  of  attack  is  a  very  practical  control 
variable  since  it  can  be  measured  and  controlled  fairly 

4  *  '' 

easily,  and  its  derivatives  are  not  present  in  the  state  or 
adjoint  state  equations.  A  control  constraint  in  the  form 
of  a  penalty  function  is  used  in  this  presentation*  This 
constraint  insures  that  the  pitch  rate  of  the  vehicle  re¬ 
mains  small  by  limiting  the  angle  of  attack.  Also,  limit¬ 
ing  the  angle  of  attack  insures  that  the  vehicle  operates 
in  the  linear  regions  of  the  lift  and  drag  coefficient 
curves.. 

A  listing'  of  the  digital  computer  program,  including 
appropriate  comments,  is  included  in  this  paper.  To  avoid 
storing  the  atmospheric  and  aerodynamic  data  on  tapes,, 
polynomial  least  squares  curve  fits  of  the  data  are  used. 
This  definitely  shortens  the  amount  of  computer  time 
required  for  each  iteration  and  reduces  the  amount  of 
necessary  storage.  The  fourth  order  Runge-Kutta  formula  is 
used  to  integrate  the  state  and  adjoint  state  equations, 
while  the  expanded  Simpson’s  formula  is  used  to  integrate 
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the  penalty  function. 

Cf  interest  in  all  such  missile  problems  is  the  evalu¬ 
ation  of  the  adjoint  state  equations,  which  requires  find¬ 
ing  the  partial  derivatives  of  lift  and  drag  with  respect 
,to  the  .state  variables,  velocity  and  altitude.  In  deriving 
these  derivatives  it  is  assumed  that  (,1)  the  vehicle  is  op¬ 
erating.  in  the  linear  regions  of  the  lift  and  .drag  coeffi¬ 
cient  curves,  and  (2)  the  parasite  drag  coefficient  is  inde¬ 
pendent  of 'atmospheric  density.  Lagrangian  five-,point  dif- 
ferentiation  formulas  are  used  to  estimate  the  derivatives 
of  atmospheric  and  aerodynamic  data. 

Chapter  II  treats  the  non-linear  equations  of  motion 
of  the  missile.  Chapter  III  then  formulates  the  optimal 
control  equations  using  the  plant  equations  of  Chapter  II 
and  Pon.tryagin's  Maximum  Principle.  The  aerodynamic  deriv¬ 
atives  used  in  the  adjoint  equations  are  derived  in  Chapter’ 
IV,  which  is  included  especially  for  control  engineers. 

The1  conjugate  gradient  technique  is  presented  in  Chapter  V. 
The  results,  conclusions,  and  recommendations  are.  presented 
in  Chapters  VI,  VII,  and  VIII  respectively. 
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II.  'Equations  of  Motion 

In.  forming  the  equations  of.  motion,  the  missile  is 

treated  as  a  variable  point  mass  acted  upon  by  thrust, 

gravity,  lift,  and  drag.  Since  the  range  of  the  missile 
Is  comparatively  short,  the  earth  is  considered  flat,  and 
the  effects  of  the  earth’s  rotational  rate  are  neglected. 

The  reference  coordinate 
system  has  its  x-axis 
along  the  surface  of  the. 
flat  earth,  its  y-axis 

vertical,  and  its  origin 

below  the  point  of  initial 

thrust.  Fig.  1  depicts 
the  free  body  diagram,  of 
the  missile.  The  fixed 
x^y  coordinate  system  is  the  reference  coordinate  system 
translated  only,  such  that  the  missile  is  located  at  the 
origin;  The  motion  measured  in  the  x-y  coordinate  'system 
Is  the  ;same  as  the  motion  measured  in  the  reference  system. 
The  equations  of  .motion  are: 

jnV  ~  P  cost*-')  -  P  -  iTig  sitftf) 

,<nvir  =  l  -  n  3  wit)  +  f  sin(oc) 
fi «  v  $y(v  (1) 

P  r:  V  C0s{t) 

4 


4 


Fig.  1.  Free  body  Diagram 
of  Missile 


GGC/EE/70-5 


where 


Appendix  B  contains  a  discussion  and  analysis  of  -Jqs  (2) 


and  (3). 


The  state  variables  are  defined  as 


where  all  the  variables  are  functions  of  time. 


(5) 


Thrust  and  Mass 

Thrust  can  be  written  as  a  function  of  time^-  Fig.  2 

is  a  graph  of  thrust  versus  time.  The  specific  impulse 

I  (Ref  10:  152)  is  defined,  in  Eq  (6)  as 

•r  _  _ .  Thrust  (F)  _  / ' \ 

Asp  “  weight  rate  flow  of  the  propel! anFs^wT  v  ‘ 


where 


w  =  %  5 


(7) 


UsSng  Eqs  (6)  and  (7) ,  the  mass  of  the  missile  can  be 


written  as 


m 


t  F 


XT* 


Jt 


(8) 
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whore  M  is  the  initial  value  of  m  —  specifically 
o 


(9) 


v?here  W  is  the  initial  weight  of  the  rocket.  Pig.  ,3  is  a 
o 


graph  of  m  versus  time.  For  this  particular  missile  WQ 
is  3000  (lbs)  and  I-  is  240  (sec). 

Equations  of  thrust  and  mass  as  functions  of  time  are 


located  in  the  computer  program  (Appendix  A)  c. 
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III.  /Formulation,  of  the  Optimal  Control  Problem, 

Cost  Function 

The  first  consideration  is  the  cost  or  objective 
function  which;  must  be  realized  before  the  optimal  control 
problem  can  foe  set  up.  The  most  obvious  cost  functions  are 

’  * 

and  (10) 

'  J= 

since  the  objective  of  this  investigation  is  to  maximize 
range.  However,  as  will  be  discussed  in  Chapter  V,  some 
penalty  on  control  is  heeded.  Therefore 

J=  ~  & 

0 

and  (11)? 

j  -  &!?$($)]■  **  ^  2  ^ 

are  used. 

Maximum  Principle  Necessary  Conditions 

It  is  desired  tp  maximize  J  subject  to 

i  -  f  ( % oL,  m,  F,  D,  L,f,t)  (12) ) 

where  Eq  (12)  represents  Eq  (1).  For  this  problem  the 
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necessary  conditions  of  Poritr, vagin*  s  Maximum  Principle  are 

(0  t  =  jf 

(*)  i  =  -  Hr 


where 


(V  nu  »  o 

&  Ti_  —  (fixfJ  fW/il  "h/ot) 

)y.  tz-bf. 

f)  tr  Jr  —j-j^F  COvfa)  S^^)j 

+  +  F?iwW  -wj  +4EM^ 

+  4  pf,  c<w  &j] 

r  -  4  IP  -  iMJ  _  F.  f,v<K)  + 

4<  -  W  J-X,  i*  [  JX,  X*  ' 


•£j  SiV(K)  + 


%  CK  fcjl,  -  *»  "  % CW 

X/  ** 

4  =  4  j  cos  (Xj,)  -  4  ■£-  - 4  4  ■<w^z) 


_  3,  21  _  2a.il 

"s  "  w  IXj  mX,W 3 

j,  =  « 

Uu  =  -«*  -  a,f|  ww  +  ¥ 
+  4E#X(  «5  W  +  g;  V] 
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Using  necessary  conditions  (1),  (2)»  and.  (3)*'  can  be 
shovm  that 

JJL  -  ft  COA/'st'M't"  '(130. 

where  the  Hamiltonian  is  not  a  function  of  time  on  the 
optimal  trajectory.  (Note:,  in  this  problem  the  Hamiltonian 
is  not  a  function  of  time  where  thrust  is  constant.) 

Boundary  Conditions 

The  specified  boundary  conditions  are 

x,(p)  -  ciin  #•/*«• 

Xi(°)  ~  0 

X-i  (0)  =  7.0,  000  -ft-  W:) 

Xs(h)  ~  0 

%  (0)  ~  O’ 

All  other  final  conditions  including  final  time  are 
unspecified*  To  avoid  adding  another  penalty  term  to  the 
cost  function  the  only  end  condition  x^,(ty )  =  0  is  relaxeu* 
Necessary  condition  (4)  is  used  to  derive  the  end 


(15) 

Of  Xtfo)  ' 
n 


conditions  .of  the  costates:, 

-  =  o 

m  =  ° 
m 6  0 

hfy  =  > 

»■ 

n  . 


"  5v.v7V  c 


./*•  VJ>  -v*’v  v*-*u*t  lyi  ■  4?  "« J  W 1 
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iV. ’  Aerodynamic  Derivatives 


In  order  to  solve  the  costate  equations,  it  Is 
necessary  to  find  the  partial  derivatives  of  lift  and  drag 
vr-lth  respect  to  velocity  and  altitude.  (In  the  equation 
lift  is  divided  by  velocity.)  Using.  Eqs  (.2)  and  (3) 


c%f'/Xi) 

)h 
LL 


~  <1*  [Ci 


ex. 


U  x*  *  v 


•i 


Mil 
XX,  J 


=  ?s 


I  ix 
z  >4 


[it  +  <  H, 


•]  *  (  x,  S  Cp 


IL  - 


•where  C. 


- 

^3 

arid  C 


!lcil  H. 

X,  f  >  Xj 
*  Cp0 


4- 


- 


+•  ^ 


2 


<s$  -i  3  Cl* 

X,  X  Xj 

Jft 


(16) 


^•]  +  ^*5Cp 


X?  f  Aj 

are  functions  of  Mach  number.  Since  Mach 


•L*  Do 

number  is  a  function  of  velocity  and  altitude  *  and  atmos 
pheric  density  is  a  function  of  altitude,  (Ref  hi,  .477') 

)  £p„  }  (p„  X/W 

77T=  x/n  xx. 


Xfp,  _ 

2^ . 

jf/w 

)X3  ■ 

XX j 

>4* 

L  1A 

(17) 

XX,  ■ 

■  XX, 

Ku 

.  M 

XX*  ~ 
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where. 


3  Cpg 

)  e 


0 


(18) 


This  assumption  implies  that  parasite  drag  is  directly  pro¬ 
portional'  to  atmospheric  density  if  velocity  is  constant. 

_ Wind_tinin_eJL_a-ri.ci Xligh t  t_e st_s_h.ay e  shown  that  this  assump¬ 
tion  is-  valid. 

By  definition  Mach  number  is  the  ratio  of  the  true  air 
speed  to  the  , speed  of  sound,  therefore 

M  i 


is  _  m  iV*  _  Vk- 

i  *3  =  Vs  )Xi  ~  Vs*  > 

where  V_  is:  a  function  of  altitude.  Substituting  Eq  (.19) 

D  '  *  ' 

into  Eq  (17), 


)  % 

m 

- 

iX  -i  ~ 


J 

‘  Vs  im 

=  iWk. 
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Curve  Fitting  Data 

A  set  of  equations  approximating  atmospheric  density 

(Ref  2:  15)  are  available  in  the'  computer  program 

(Appendix  A).  Fig.  4  is  a  graph  of  density  versus  altitude. 

VB.  CDo-  and  C^  were  approximated  by-  using  the  piecewise, 

*  < 
polynomial,  least-squares,  curve-fit  method  (Appendix  D) . 

Tne  curve  fits  are  graphed  in  Figs.  5t  6,  and  ?.  The 
curves  with  the  asterisk  are  the  curve  fits;  the  curves 

without  the  asterisk  represent  the  given  data  points. 

Data  for  the  speed  of  sound  were  .obtained  from  Ref  (8:  4). 

Cr,  and  Cr  ,  data  were  obtained  from  wind  tunnel  tests. 

Do  L- 

Table  I  shows  the  maximum  deviation  in  percent  between  the^ 
ordinate  of  the  given  data  points  and  those  of  the  curve 
fit.  The  curve  fit  polynomial  equations  are  also  located 
in  the  computer  program  (Appendix  A). 

Table  I 

**  Errors  in  Curve  Fits 


Dependent 

Function 


Maximum  Deviation 
in  Per  Cent 


Value  of  the 
Abscissa 

3.5  x  10^  (ft) 

3  (Mach  No.) 

1  (Mach  No.). 


^Differentiating  the  Data  Curves 


In  order  to  solve  Eq  (20)  and  then  Eq  (16),  it  is 
necessary  to  find  the  partial  derivatives  of  V_  and  f  with 

S 

respect  to  altitude,  and  those  of  CLet  and  Cp-0  with  respect 
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to  Mach  number.  An  approximation  of  these  derivatives  is. 
obtained  by  ysing  the  average  of  a  set  of  five,  five-point, 
Lagrangian  differentiation  formulas.  A  derivation  of  the 
formulas  is  located  in  Appendix  C.  The  approximations  for 
the  errors  in  the  formulas  are  assumed  negligible. 

It  is  now  possible  to  solve  Eq.  (1.6)  if  the  velocity 
and  altitude  of  the  vehicle  are  known. 


c 
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Fig*  Atmospheric  Density 
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V..  The  Conjugate  Gradient  Technique 

Nov;  that  the  optimal  control  pro  olem  has-  been 
formulated,  some  iterative  technique  must  be  used  to  solve 
it.  In  this  paper  the  conjugate  gradient  method  ‘(Ref  7),  is 
used.  In  general  this  method  is  able  to  converge  quickly? 
to  a  near  optimal  solution  from  ;poor  initial  guesses  on 
control . 

The  Conjugate  Gradient  Algorithm 

(1)  Choose  and  arbitrary 

(2) .  Set  pQ  =  Hu(o*0)  where  i  =  0 

>•(3';)’  Find  k.^  such  that  +  k^p^)  is  maximized 

with  respect  to 

(4)  Set  <p(iFl  =  +  k^/ 

(5)  Set  Hui+1  =  Hu(^i+1) 

(.6)  Test  the  gradient  for  convergence;  if  there  Is 
convergence,  stop. 

(7)  Set  ^  =  ^Hui+1*.  Hu^.  ( inner  product) 

(8) '  .Set  *P^+2.  =  ^ui+l  **  @  i^i 

(9)  Repeat  starting  with  Step  (3)* 

It  is  desired  to  f ihd  <*(i+1  such  that  Hu1+1  is  zero,  for 
all  time.  Because  this,  problem  is  non-lineai*  and  non- 
quadratic,  the  conjugate  .gradient  method  does  not  find  the 
optimal  control  in  five  Iterations.  Therefore,  it  is 
necessary  to  set  some  arbitrary  tolerance  oh  convergence 
and/or  stop  the  program  using  some  other  condition  such  as 
the  number  of  iterations  or'  computer  execute  time.  (The 
gradient  in  Step  (7)  is  treated  as  an  n  by  1  vector  where  . 
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M  =  p  +  1-). 

The  Linear  K~ Sear ch" Algor i thm 

The  linear  K- search  j.s  the  method  used  to  perform 
Step  (3)  of  the  conjugate  gradient  algorithm.  It  is 
assumed  that  J  is  a  linear  function  of  cl  and  therefore, 

=  p^Hu(c<^+lCjP^)’ .  In  this  problem  J  is  definitely  ■ 
not  a  linear  function  of  but  it  is  assumed  that  ,p^Hu  is 

a  good  approximation  of  j( 

* {  T K-  <-  )  “  TfrU)} 

consider  p^  and  Hu^  n  by  1  vectors. 


(1)  Initial  guess  of  k, : 


(i)  '  Let  K  =  iK/if  (KK^Pj^r1 

(ii)  Let  K  =  if  is  elsewhere 

(2)  Evaluate  J (o^-hnKp^ )  where  (n=0,l,2,4, 8, . . .  a.b) 
and  PiTHu(«<^4-bKpi)-«CQ;  therefore  aK^k^bK 

(3)  Interpolating  k^: 


(i)  Let 


where 


“■  Jf/(1  irr1-  i-^+Ytu 

W  =  fe*-r  Jtaj  7 Ik)]  2 

J(a)  r  7  (°*  a  +  ft  K  Pj )' 

J  f&J  z-  fj ik  'fai + 4- 
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(ii)  Set  ^1+1  =a'i  +  k^r,  if  J(a)  and 

J(b),<  JfkjO 

(iii )  If  neither  d(a)  nor  .:{?)  is  less  than 
J(k^),  return  to  See;:  .2)  using  n  =  a  + 

where  (  j  =  1, . . .  =(  .  .5)  and 

Jk<  kJL<  ik 

With  the  exception  of  Step  (3-iii this  linear  search 
is  the  linear  search  method  of  Fletcher  and  Reeves-  (Ref  3)  • 

The  equivalence  to  Step  ('3-iii)  in  Flerfner  and  Reeves 
divides  the  interval  into  two  subinterreis  at  k^  and  tests 
the  sign  of  3"(k^);  to  determine  which  sv.: interval  to  use  in 
Step  (3-i).  The  Fletcher  and  Reeves  me -.hod  is  less  cumber¬ 
some,  but  Step  (3-iii)  seemed  more  reliable  in  this 
problem. 
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VI .  Results-  and  Discussions 

The  maximum  range  obtained  in  this  investigation  is 
sensitive  to  the  guesses  of  the  fixed  final  time,  the 
weighting  function  R,  and  initial  control.  The  .-selection 
of  values  for  these  three  parameters  determines  whether  a 
local  or  global  maximum  range  is  obtained.  A  trial-and-* 
error  method  is  perhaps  the  only  way  of  etermining 
reasonable  estimates  for  the  optimal  values  of  those 
parameters  *  - 

The  Final  Time 

Some  idea  of  the<  values  for  the  fixed  final  time  was 

obtained  by  computing  trajectories  with  constant  angles  of 

attack.  The  values  used  were  0.05*  0.1,  O.15,  and  0.2 
radian..  Judging  from  their  impact  time,  it  was  decided  to 
work  with  tf  =■  300  seconds. 

The  Weighting  Function 

l'n  order  to  get  some  idea  of  the  best  magnitude  of  the 
weighting  function,,  conjugate  gradient  computer  runs  were 
made  with  R  and  equal  to  zero  for  all  time,  and  tf  equal 
to  various  values  above  300  seconds,  in  general  these  runs 

were  unstable,  i.e.,,  on  many  iterations  control  would 

greatly  exceed  20°,  the  flight  path  angle  would  continuously 

increase,  and/or  the  range  would  be  negative.  However,  on 

the  more  stable  runs  the  values  of  the  gradient  gave  an 

_  Ji.  _  _  _ _ 23.  .  . 
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indication  for  the  necessary  magnitude  of  R.  For  cost 

h  < 

c^e  gradient  had  magnitudes  around  10  and  lO-' 

'(lO^  and  during  thrusting,  while  during  glide  the 

magnitudes  were  around  10^(10^)  and  lower.  Also,  the  angle 
of  attack  would  exceed  20°  during  thrusting*  specifically 
during  the  first  10  to  15  seconds,,  and  as  expected  would 
remain  well  below  20°  during  glide.  In  order  for  the 
penalty  function  to  have  some  .affect  on  the  gradient  during 
thrust,  the  magnitude  of  R  would  have  to  be  comparable  to 

that  of  the  gradient.  The  function  used  in  this  presenta¬ 

tion  was 

r  =  (l0  [  ait)  -  n(* -*<■)];  <21> 

k  £ 

where  RQ  is  some  number  in  the  neighborhood  of  10  and  10^, 

u(t)  is  the  unit  ste,p  function  and  tr  is  the  approximate 

complete  burnout  time  of  30.00  seconds.  (In  the  following 

discussions  t  will  take  on  different  values.) 
r 

The  Initial  Guess  on  Control 

Once  some  idea  of  the  magnitudes  of  R  and  t^  had  been 
obtained,  the  next  step  involved  checking  the  sensitivity 
of  the  conjugate  gradient  method  to j various  initial  guesses 
on  control  Four  guesses  (0.0,  0.1,*  0.15,  and  0.2  radian 

for  all  time)  were  tested  by  using  and  varying  RQ  and' 

until  what  appeared  to  be-  the  best  sample  runs  for  each  of 
the  four  values  was  obtained.  On  the  basis  of  the  sample 
runs  =  0.15  radian  appeared  to  be  the  best  choice. 

l- 


4' 
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However,  each  of  the  four  'guesses  had  one  major  fault:  the--' 
conjugate  gradient  method  did  very  little  in  optimizing 
control  as  time  approached  t^ ,  consequently,  at  final  time 
large  negative  flight  path  angles  were  obtained.  Therefox-e, 
the  initial  guess  on  control  was  chosen  to  be 


£>.}£;  tC 

0.00l(t~W)  +  0>l$J 


t  >  3i& 


(22) 


was  increased  linearly  as  time  approached  tf  so  that 
y(tf)'  would  be  increased,  thus  increasing  range.  The  time 
of  315  seconds  was  chosen  in  Eq  (2£)  because  in  that 
neighborhood  the  flight  path  angle  reached  critical  nega¬ 
tive  values.  It  must  be  mentioned  that  many  other  functions 
for  UQ  where  t  >  315  seconds  could  have,  been  used,  and  even 
the  use  of  t  -  315  seconds  as  the  break  point  is  question¬ 
able. 


■  -Results  of  Sample  Computer  Runs 

Using  Eqs  (22.) ,  (2/),  and  the  intervals  1  x  10^  ^* 

6  c 

,R0  1\  x  10  and  300  ~  t^  ^  380<  sample  runs  were  made  to 

determine  smaller  intervals  for  R„.  and  t».  If  the  run  were 
-  -  or 

unstable,  then  R„  was  increased  by  an  amount  that  was  some- 
what  proportional  to  the,  amount  of  instability.  If  the 
altitude,  at  time  t^.  were  large,  then  t^  was  increased  so  as 
to  decrease  h(tj.,),.  Table  II  on  the  next  page  depicts  the 
results  of  some  of  the  more  significant  runs.  In  the  table. 

.  25 


GGC/EE/70-5 


a  sample-  run  was  one  with  a  computer-execute-time  limit  of 
four  or  five  minutes*  These,  runs  were  used  to  determine  if 
a  complete  run  using  the  sampled ‘guesses  was  warranted.  A 
complete  run  was  forty-one  conjugate  gradient  iterations. 

Using  the  above  procedure  the  intervals  were  narrowed 
to  4.5  x  10^  £  R  ~  9  x  10**  and  360  t^.  &  380.  Runs  (1) 

through  (8)  of  Table  II  depict  the  results  obtained  using 

R  and  t~  within,  those  intervals.  Run  (1)  obtained  the 
o  1  . 

✓ 

maximum  range,  however,  it  was  unstable  and  the  gradient  of 
the  maximum-range  trajectory  was  very  large.  Runs  (2.) 
through  (8)  were  attempts  to  stabilize  the  method  and  hope¬ 
fully  increase  the  maximum  range  obtained.  In  these  runs 
stability  was  obtained,  however,  the;  range  obtained  in  Run 
(l5)  was  never  equaled,  or  improved  and  very  little  improve¬ 
ment  was  made  on  the  gradient.  Also,  varying  tf  did  very 
Tittle  to-  Improve  range  maximization..  In  fact  Run  (3)  with 
t^  =  36O  seconds  was.  never  Improved  as  tf  was  increased. 
Therefore  all  succeeding  runs  were  made  with  t^.  =  3&0 
seconds. 

The  next  step  involved  decreasing  tr  over  the  interval 

5.0  «£•  tr  ^  30.  while  varying  RQ  over  the  interval  1  x  10-* 

Rq  —  8  x  10^.  The  best  results  were;  Runs  (9)  and  (10=) 
where,  the  maximum;  range  obtained  was  143. 76  miles.  Still 
the  gradient  was.  large  during  thrust.  ,At  this  point  in  the 
investigation  it  was  discovered  that  the  mass'  of  the  Vehicle 
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was  being  treated  as  a  function  of  altitude,  i.e..,  in  the 
computer  program  "g"  of  Eq  (4)  was  also  being  used'  as  Mg  ,- 
of  Eq  (9) i  Although  this  error  affected  the  range  of  the 
trajectory  by  tenths  of  a  mile,  the  value  of  the>  best  Rq  was 
greatly  affected.  A  comparison-  of  Runs  (9)  and  (11)  gives 
an  indication  of  the  affects  of  the  correction.  Runs  (Hi) 
through  (.14)  used  the  correct  value  for  ‘"jg  P, 

Using. -a  modified  history  of  the  control  obtained  in  Run 
(9)  as  the  initial  guess  (.this  control  affected  a  range  of 
143 ;26  miles,-  t  =  3&0  sec. )  various  runs  were  made  in  an 
attempt  to  improve  the.  range.  Runs  (12)  through  (14)  are 
samples  of  those  runs.  Because  of  the  results  of  Runs  (2) 
through  (8),  no  complete  runs  were  made  and  t  =  10  seconds 
was  used.  As  can  be  seen  very  little  improvement  was  made 
on  the  maximum  range.  =Further,  the  gradient  was  not 
improved . 

At  this  point  it  was  decided  to  use  and  the  control 
history  obtained  in  Run  (14)  as  the  initial  guess  on  con¬ 
trol.  A  graph  of  this  control  history  is  depicted  in  Fig. 

8.  Using  tf  =  3'60  seconds,  R  was  varied  from  0.1  X  10 
14 

to  0.1  X  10  using  the  same  scheme  -that  was  used  with  J^. 

No  improvement  on  the  range  was  obtained.  Most  of  the  runs 
were  unstable;  however,  the  more  stable  runs  were  unable  to 
search  effectively  in  the  increasing  range  direction.  Con¬ 
sequently.,  the  conjugate  gradient  approach  was  abandoned 
even  though  the  optimal  trajectory  had  not  been  obtained. 


Control  His 


GGC/EE/?0~5 


Two  Modifications 

Using  the  control  h,istory  of  Run  (14),  two  modifica¬ 
tions  were  made  (1)  different  constant  control  values 
were  simulated  for  the  time,  region  t  >  65  seconds,  and  (2) 
the  conventional  technique  of  flying  the  minimum  drag-to- 
lift-ratio  trajectory  during  glide  was  simulated.-  As  a 
result  the  trajectory  was  improved.  The  control  history  of 
Run  (14)  obtained  a  total  range  of  1*17.64  miles  at  t  = 
385.73  seconds  (control  was  held  constant  at  11.46°  after 

V 

t  ~  360  seconds).  In.  modification  oner  e<  -  0.2,  ,0.225, 

0.25,  and  0.275  radian  were  used.  tk.  -  0.25  obtained  the 

best  range  — •  156.79  miles  at  t  =  450  seconds. 

*  '  ’ 

The  necessary  condition  for  modification  two  is  that 

i  iv 0)  _  Q 

)  *  (23) 


Substituting  Eqs  (2)  and  (31  into  Eq  (23), 


(24) 


where satisfies  Eq  (23)?.  Using  the  trajectory  obtained 
in  modification  one  and  o(  =  two  trajectories  were 

simulated  —  the  first  used  U  =  LD  for.  all  of  the  glide, 
and.  the  second  used  it  for  t  >  65  seconds.  The  best  range 
obtained  was  154.27  miles  where  o(  =<?>££)  f ° t  >  65  seconds,. 


The  Best  Trajectory  Obtained 

The  best,  control  history  (Fig.  9.)  in  this  investiga¬ 
tion  obtained  a  range  of  15.6.79  miles.  This  history  is 
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well  within  the  limit  of  20  since  its  maximum  point  is 
16.91°.  The  control  programs  of  runs  in  Table  II  that 
reached  at  least  140  miles  had  the  same  general  shape.  Those 
trajectories  that  fell  below  140  miles  had  control  histories 
with  larger  minimum  values;  consequently,  higher  altitudes 
were  obtained.  As  expected  the  range  is  very  sensitive  to 
the  angle  of  attadc  during  thrust  and  slightly  thereafter. 

Figs.  10,  11,,  12,  13  and  34  are  graphs  of  the  -state 

variables.  Fig.  14  includes  a  graph  of  the  pitch  angle  0. 

Fig..  15  is  a  graph  of  the  pitch  rate,.  The  maximum  acceler- 

2 

ation,  4.37gQ»  end  the  maximum  pitch  rate,  0.047  rad/sec 
occur  around  20  seconds.  Both  values  are  quite  coierabl'e 
for  equipment  design  purposes.  The  pitch  angle  at  final 
time  Is  equal  to  -2.64°  ahd  remains  within  -3*7°  and  -2.2 5° 
for  the  last  50  seconds  of  flight.  This  is  quite  adequate 
for  llne-of- sight  and  line-of-slght  rate  steering.  The 
velocity  at  impact  is  Mach  0.6. 

Figs,.  16  and  17  are  graphs  of  the  gradient  and  the  Ham¬ 
iltonian.  Both  graphs  indicate  that  the  optimal  trajectory 
has  not  been  obtained..  Note  that  -as  -time  approaches  final 
time  the  gradient  is  small  and  the  Hamiltonian  is  almost 
constant.  These  facts  indicate  the  reason  that  the  conju¬ 
gate  gradient  method  is  inactive  in  this  region,.  Also  evi¬ 
dent  is  the  influence  of  the  interval,  19*44  £  t  30»35t 
where  the  Hamiltonian  is  a  function  of  time.  In  these 
graphs  R  =  0.  • 
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VII.  Conclusions 

A  cost  function  including  an  integral  penalty  function 
to  constrain  control  can  be  used  to  obtain  an  approximation 
cf  the  range-maximizing,  angle-of-attack  control  history  of 
an  air~tc- surface  missile  trajectory  using  non-linear  equa¬ 
tions  of  motion  and  the  conjugate  gradient  method. 

The  major  shortcoming  is  the  necessity  of  having  to 
guess  values  for  the  weighting  on  <the  penalty  function,  the 
initial  guess  on  control,  and  the  final  time.  The  maximum 
range  obtained  is  greatly  dependent  upon  these  guesses. 
Another  shortcoming  is  the  inability  of  the  method  to  opti¬ 
mize  control  as  t-ime  approaches  the  guessed  final  time. 
Finally,  the  most  obvious  drawback  is  that  the  exact  opti¬ 
mal  trajectory  is  not  obtained. 
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VIII.  Recommendations 

Although  the  approach  used  in  this  investigation  is  one 
of  the  simplest  ways  of  setting  up  the  problem,  obviously', 
other  ways  might  lead  to  better  results.  With  the  exception 
of  the  last  two  subsections,  all  of  the  following  recommen¬ 
dations  discuss  possible  ways  of  improving.  \the  results  of 
this  paper*.  The  subsection  "Numerical  Methods"  treats  pos¬ 
sible  ways  of  decreasing  the  computer  execute  time  necessary 
for  each  gradient  iteration,  and  the  last  subsection  "The 
Control  System"  -treats  a  possible  area  of  further'  investiga¬ 
tion 

Objective 

As  with  all  optimal  control  approaches  there,  are  a 

i 

number  of  ways  of  setting  up  the  problem.  ’  As.  suggested  by 

i 

Las  don-;  Matter,  and  Waren  (Ref  7)  end  conditions  can  be 
treated  using  penalty  functions  in  the  objective.  Since  it, 

is.  desired  that  x^(t^).  =  Oj  the  objective  can  be'  written  as 

^  + 

7-  -K,  [XM)]^  [W]‘-  A)a  yx^  (26) 

where  R^  is  a.  weighting  function.  Thus,  the  use  of  this 
objective  tends  to  maximize  the  range  while  minimizing  the 
altitude  at  time  tf.  However,  this  approach  is  very  sensi¬ 
tive  to  the  values  used  for  the  final  time  tf  and  R^.  Of 
course  trial-and-error  methods  would  have  to,  be  used  to  de¬ 
termine  reasonably  good  values  for-  tf ,  R^,  abd  R. 
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A  better  objective  would  be  ^ 

t  ■ «,  fabft  [y,hf-  4^Ji: 

which  tends  to  maximize  range  as  well  as  altitude  at  final 
time.  This  approach  is  less  sensitive  to  final  time. 

Constraining  Control 

The  problem  of  constraining  the  control  can  be  handled 
by  introducing  a  fifth-,  state  and  dropping  the  penalty  inte¬ 
gral  term-  of  'Eq  (10).  Since  it  is  desired  that  4» 

,  then  the  fifth  state  equation  can  be 

Xf  -  (28O 

This  keeps  the  objective  as  simple  as  possible,  but  adds 
another  state  variable  and  adjoint  state  variable  to  the 
problem.  Again,  "brute  force’1'  techniques  will  have-  to  be 
used  to  determine  R^.. 

Method  of  Second  Variation 

A  more  sophisticated  optimal  computational  scheme  is 
the  second  variation  method,  (Ref  9:  4l4).  The  conjugate 
gradient  method  essentially  searches  for  the  first  order 
effects  of  the  control  on  the  objective.  By  considering 
second  order  .effects;  as  well,  the  second  variation  method 
converges  much  more  rapidly  than  the  conjugate  gradient 
method  (Ref  7:  138).  However,  the  second  variation  algo¬ 
rithm  is  much  more  complex  and  very  sensitive-  to  the  initial 
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guess  on  control.  (In  general  if  the  initial  guess  on  con- 
trol  is  not  close  to  the  optimal,  the  second  variation 
method  will  diverge . )  Of  interest  is.  the  use  of  the  angle- 
of-attack  history  obtained  by  this  investigation  as  the 
initial  .guess  on  control  in  the  second  variation  method. 

Two  Control  Variables 

Another  sophisticated  approach  is  the  use  of  the  conju¬ 
gate  gradient  method  in  solving  an  optimal  control  problem 
involving  'two  control  variables,  specifically  the  angle  of 
attack  and  the  thrust-vector  angle.  (The  thrust- vector*' 
angle  is  measured  from,  the  x  -axis  of  the  missile  to  the 
thru  at  vector.  In  this  paper  the  thrust-vector  angle  is 
zero  for  all  time.)  During  thrusting,  the  gradient  becomes 
a  2  X  1  vector  and  the  conjugate  gradient  k-parameter  be¬ 
comes  a  2  X  2  matrix.  After  -thrusting,,  the  problem  simpli¬ 
fies  to  a  one-control -variable  problem.  Mo.st^if  not  allj 
the  matrix  subroutines  needed  to  program  this  two-control- 
variable  approach  are  stored  on  the  7044/7094  computer 
system  library  of  the  Digital  Computation  Center;  WPAFB, 
Ohio..  ' 


Interval  Maximization 

One  major  fault  in  the  approach  used  in  this  investiga¬ 
tion  is  the  tendency  for  the  conjugate  gradient  method  to 
make  little  or  no  effort  to  iterate  toward  an  optimal  tra¬ 
jectory  in  the  region  where  time  approaches  the  final  time. 
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Maximizing  over  smaller  arid  smaller  intervals  of  time  seems 
to  be  a  way  of  solving  tfcls  problem.  First  maximization  is 
attempted  over  the  entire  range  of  time  (t0,  tf)  as  was  done 
in  this  paper.  Then  some  time  t^  is  chosen  where  t-  repre¬ 
sents  the  time  beyond  which  the  conjugate  gradient  method  ap¬ 
peared  to  fail.  Next  maximization  is  done  over  the  interval 
(t^»  t-.) ;  Then  a  larger  is  chosen  and  maximization  is 
repeated  over  that  smaller  tvime  interval. 

Numerical  Methods 

As  previously  mentioned,  a  set  of  Lsgrangian  differen¬ 
tiation  formulas  is  used  to  calculate  the  approximations  of 
the  derivatives  of  the  data.  Using  these  formulas  in  the 
computer  program  is  somewhat  cumbersome  arid  time  consuming. 

A  more  efficient  way  is  to  use  curve  fits  of  the  approxima¬ 
tions  of  the  derivatives.  However,  curve  fitting  the.  derdv- 
'atives  is  a  more  delicate  operation  than  curve,  fitting  the 
given  data  points.  (An  effort  to  curve  fit  the  approxima¬ 
tions  of  the1  derivatives  with  polynomials  resulted  in  curve 
fits  that  are  not  as  accurate  as  those  obtained  in  Chapter 
III  and  consequently,  are  not  used  in  this  presentation.,) 
Judicious  use  of  a  combination  of  least  square  curve  fits 
such  as  the  polynomial:  arid  the  exponential  may  result  in 
curve  .fits  with  satisfactory  accuracy. 

In  this  presentation  the  fourth  order  Runge*?Kutta  inte¬ 
gration  formula  Eq  (1A)  is  used' to  integrate  the  entire 

missile  trajectory.  .Although  very  accurate  and  quite  stable, 
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this  method  is  somewhat  cumbersome  and  time  consuming. 

Other  methods,  such  as  the  Adams-Bashforth  predictor-correc¬ 
tor  method,  can  be  just  as  accurate,  but  far  less  time  con¬ 
suming,.  The  Runge-Xutta  formula  can  be  used  to  determine 
the  first  four  integration  points  and  then  a  faster  method 
can  be  used  to  determine  the  succeeding  .points*  Of  course  a 
more  sophisticated  integrating  subprogram  is  necessary,  but 
computer  execute  time  is  reduced. 

The  Control  System 

Once  the  optimal  sngle-of-attack  history  has  been  ob¬ 
tained  a  further  investigation  treats  the  design  of  a  prac¬ 
tical  control  system  that  flys  the  missile  along  the  optimal 
trajectory.  In  designing  such  a  control  system*  open  or 
closed  loop,  some  of  the  aspects  that  must  be  considered  are 
the  desired  accuracy ,  the  weight  of  the  control  system,  the 
space  available  on  board  the  missile,  the  cost,  the  type  of 
control  -system*,  the  possible.,  use  of  other  control  variables 
besides  angle  of  attack,  and  the  fact  that  the  system  will 
have  to  operate  in  real  father  than  standard  atmospheric 
conditions.  Finally,  one  decision  that  must  be  made  is 
whether  to  use  classical,  optimal,  or  stochastic  control  de¬ 
sign  techniques.  All  three;  have  their  own  advantages  and 
disadvantages . 
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Appendix  A 


Computer  Program 


The  following  computer  program,  written  in  Fortran  IV, 
takes  on  the  average  less  than  0.5  minutes  of  execute  time 
Ctf  =  3-60 ;  t  =  5)  to  perform  a  conjugate  gradient,  iteration 
on  the  IBM  70*1-4/7094  II  Direct  Coupled  Operating  System  of 
the  Digital  Computation  Division:  Aeronautical  Systems  Divi¬ 
sion,  ,W  r  i  gh  t  -  Pat  t  e  r  s  on  AFB,  Ohio.  The  program  is  composed. 


of  twelve  subprograms: 


(1) 

MAIN 

*■* 

performs  the  conjugate  gradient  algo¬ 
rithm 

(2) 

NEWCON 

- 

performs  the  k~ search,  consequently 
finding  the  new  guess  on  control' 

(3) 

GRADNT 

- 

computes  the  gradient  . 

(4) 

EQUAT 

uses  the-  Lagrangian  formulas  to  differ¬ 
entiate  the  aerodynamic  data;  contains 
the  state  and  adjoint  state  differen¬ 
tial  equations 

(5) 

PMAS 

- 

contains  the  equations  for  the  mass  of 
the  missile, 

(6) 

PTH3US 

- 

contains  the  equations  of  thrust 

(7) 

FRHO 

- 

contains  the  equations  approximating 
atmospheric  density 

(8) 

FVS 

— 

contains  the  polynomials  approximating 
the  velocity  of  sound 

(9) 

FCLA 

- 

contains  the  polynomials  approximating 

the  C,  .  curve  of  the  missile 

L  cc 

(10), 

FCDO 

- 

contains  the  polynomials^  approximating 

the  C-.  curve  of  the  missile 

Do 

(11) 

INT3G 

»• 

uses  Runge-Kutta  fourth  order  formula 

to  integrate  the  state  and  adjoint 
hh. 
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state  equations;  calculates  the 
objective 

(12);  PRICE  -  uses  the  expanded  Simpson  integration 
formula  to  integrate  the  penalty 
function 

The  Runge-Kutta  formula  used  is 

XlKl  =  X„  +-  j(  !<,  +  2  •!'  2 H  +  Z'l)  <1A) 

where 

=  dffxftn,) fn)l 

Kg  =  Xn  +  r)] 

it,  C  f,  +  &)  ] 

The  Simpson  integration  formula  is 

£%ft)  =  f(t>+ H**'i**  V  ■  ■  ■  +HH+t,0  (2A) 

' 


where 
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I 


5  I  9  JCR 
cniFT-C  MAIN 


MAf 


I  *  * 


§:• 


l-'O 


i-i 


<4  I  N 

..com  ion/too/  dt  imp  ♦  .v-  rr.r ,  nst  Air «  sarfa  » cost  *  fti  >  \ ,  fcontr  ♦  xvosr 4  «o  4  t i 
D  I  MEN?  I  ON  X  (F  4  1  00 )  I'CCMTR  ( ICO ) *  T I  MF  (  ,100  )  4  OHDU  (100)4  DHDUU  1  00  )  » 
10(100)  4A(  100)  flt  lfiO) 

READ (  5  i'l  CO )  NST!™  .  NOTATE  *  I  TEH  »DT  !f‘E  « TOt. » T 1 1 .00 
FORMAT  (OIL  i'3r  1 0  *  :■  i  F 1 0  i .?  > 

1  JX--0 

SAPEAsff. la 
NSTAt  ?..r:r  TAT*'v  /. 

NSTAT 1  -NST  AT 1 

N  ?.  T lr  c‘  3  "  Iv'  0  TrP'  -•  0 

RF  /.n  (5,101)  ( X  (1  4  1)  »  1  =  1  » NST ATF )  <  ( X  (  1  1  N5.TEP3  >  4  1  =NST AT  1  4 NSTAT2  ) 
FORMAT  Ufl'.C  ) 

7  I  MEM  )-a.o 
’NSTEP  1  “NATE^W-  1 

do  li  !-s<Nsr.r:-  i 


I  M 

F  (1 

1  L 

+  DTI  Ml 

If. 

*  6 1  / 

37 

2‘  nf< 

19 

».37/ 

57 

2938 

10 

.  37/ 

37 

2958 

7. 

9>7  / 

57 

9958 

0 

.  77/ 

87 

/L  >  r-  L 

-1 

»  99/ 

87 

1 

♦  or./ 

r>7 

r?9rC 

A 

•  65/ 

87 

1C7 


102 


CON TO <  O ) =  6 . 6 3/07 . 29  rf! 

CON  Vi-  (ICO  "7 . 62/07 «  2900 
CONTR ( 1 1 )"S*  1 0/57 «  2950 
C ONTR  (  1 2 ) =8 . 35/57 .  2  35?. 

C-ON'TR  (  1 3 )  af) ,  4 8/p 7 , 292 F 

CONTR ( 1 A ) =8.61/  J7 . 2988 

DO  30  I  =  1  *~4 100 

CONTR ( 1 ) =0.35 

CALL  l NT EG ( T! ME i X , CONTR ) 

CALL  GRADNT  (7  I M~  4  X  ,  CONTR « OHDU.) 

C0CT  =  2.,0'**(  COST  5XV0  ST ) 

COST- SORT ( COST  > /52S0 » 

XVOST  =  XV037/5?.?0 .  *  *2 
7/RITE  <6 4 107)  COST « XVOST 

FORMAT  (  1 4'H  INITIAL  COST  =  i  1  PE  15*8  rl  5H  INITIAL  XVOST- 4 1  PET  5 1 8) 
XVOST = XVOST *5280 . 

COST -*0057*5286  • 

COST  =  0 .  5 ■"•COST  ”  *2~X VCST 
CO  I  Ir-L-jNcfEPl 
OHO’Jl  Ci  )-DwDU(  I  ) 

SCI)  =rvi-?»u<  i  ) 

CALL  NEW  COM  ( T  I  M“  4  X  4  CCNT^  4  S » PHQ'J  } 

I J<~ i JX+ l 
WRITE (6 4 10? )  iJK 

FORMAT  (  1  ?!  i  1  I  TER  AT  LON  NO .  4  I  5//90H  T  I  ME 
ALTITUDE  RANG”  CONTROL 

( D~G3  T"T ) 


1  G  AL  i  i  il> 

Z~~.)  ( FT  /SEC) 

DO  2  1-1 4 NS TE0 I 
X  ( 2  4  1  )  -  X  ( 2  4  I  )  *57 . 2f*58 
X(A  «  I  )=X(d  4  I  )  /52-30. 

CONTR (  1  )  -CONTR (  I  1**57.2953 


VELOC I  TV  PATH  Al 

GRADIEN7/69H  (: 

\  I  )  \  L/C.  O  ;  / 


c 


o  no  o  o  o 
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V,:R I IF  (6*103)  TIM£(  I  )  »X(  1 «  H  4  X  ( 2  4  I  )  « X  Ct> ,  I  >i  X ( A  ,  J  >  iCONTrM  1  J  *DHDU (  I  ) 
F 0PM AT  ( f  iB  «2  4  AX  *  F9 . 2 » 4X  4  F7>  ?>  4*  *  H'O .  ?i4X4  FB  *  P.  t  A  X » F  7  »  ?  4  1  OX 4 1  PE  1;5  .  B ) 
X<?«  M=X(?4  I-)/f57*?95S 

X  ( 4 « I ) «X ( A  »  M *3230 *  • 

CONTR ( I>  =CONTR < I ) /57 . 29*V 
] F ( T IMF ( i ) . GT  «  FT l K  >  GO  TO  TO 
COiVT  I  NUc 

COST  =C0ST/“T?00  *  %#?. 

HR J-TE  ( 6  ♦  1 00  )  COST 

FORMAT  (//l  1H  TRUE  CCST“»;  F8*2) 

COST =G0ST*3280  « 

F C ONTR =■ F COMT R ?.  07 . 2  7~c. 

RANGFs?  <.0~  (  CO ST vX VO SI  5 
PANGS- SORT  M?A -vG*-*  !  /BFOO. 

VP  I  T>~  { 6  4  2  04  );  RANGE  f  I M  4  FCCNTR 


FORMAT!  Tii  RANGE- 4  F8-2.4 
1L  AT  FINAL  TI?.;r>.irQc2)‘ 
FCONT  RaFCOMTP/riV  -  29D3 
XVOST  =  XVOST/5200 .  «°S'2 
WR I TE ( 6  4  1 06 )  XVOST 
FORMAT  (  7K  XVOST  "...  T  PE  i  5 « , 
XVOST  -  XVOST  •::-n28b  *«■  v2 


12H  FINAL  T'l MZ~ 4  FS  e  2  4 32H  VALUE  CF  CCNTRC 


TEST  FOR  CONVERGENCE 


TGRAOsO.O 

DO  3  \l  =  1  « N5TEP1 

BH=ARS<DHDU(I  )  ) 

I  F  ( OH  .  GT  4  f  GR AD')  TGR AD=DH 
WR  I  TE  (  6*4  1  05  )  TG°AD  4  TCL 
FORMAT (A AH  MAX#  ABSOLUTE  VALUE  OF 
119H  THE  TOLERANCE  i S 4.1  PETS *6) 

I  F  (  1  J;<  tGT «  I  TER  )  GO  TO  4 
! F ( TOL • GE  «  TGRAD )  GO  TO  4 


*E  GRAD  I  EN-T  I  S  9  1  PE  1 5  »  8/ 


FINDING  NEW  SEARCH  DIRECTION 

DH=C.. 

DH1=C. 

DO  ii  i  =  H4rMSTEP2' 

OHaOWnUftTl  **2*0H 

DH 1  =6kOU  J  <  V)  *«*2+DM  1 

BETA-OH/OHi 

DO  8  1=1 4NSTEPI 

S  <  I)  =  DHOU  (  I  TFPETA-icS  (  i  ) 

DHOUi  (  I  )=DHr>U(  I  ) 

GO  TO  9 

STOP 

EK'D 

SUBROy'T  I  \'Z  NE  V.’CCN  (  T  I  ME  4  X » CG.NTR  4  S.4  DHD'J  ) 

COMMON/TOO/  DT  I  ME  1 NSTZ?  4  NSTATEi  S  AREA  1  COST  4  #'T  I M  4  FCON'TR  1  XVOST  4  PC'  »  T I 
D  I  MENS  i  ON  X  <  S  1  1  00  io.  CONFR  {  I  CO  )  I  ME  (1 OC  >  4  DHD'J  (  J  00 )  4  S  (  j  ^0  ) 

T)0  1-33:  L- I'.rNSTEP  V  4  *  ' 

1  ViR  I  TE  <  6  4  1 30  )  T  I ;  \E  (  I  >  1  S  (  i  ) 

’J  FORM  A  T  ( F  i  0 . 2  4  EX  4  1 P  F 1 5  «  S ) 


ESTIMATE  CF  'ORDER 


MAGN I TUDE 
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C 

r- 

V.- 

C 

C 


TCOST  =  COST 
Xi<-0*6 

N’STfcP’NWSTIrP+l  .  ; 

DTC0ST=6:*d 

DO  i  1=1  4NSTFPI  "  -  -•* 

DTCOSTaDTCOST-kS  (.  I  )  *DHDU  < -I  ). 

!  >;ksxic+s(t)-hp 

XK=  1R0/SQPT(  XK) 

DO  2  !  -l  4XSTFPI 

2  CON-TPi  V)  =COMTpn  >  *  >.•  K-”-S  ("I  ) 

CALL.  I NTEG ( T  WIE-j&C-t  CONT‘5  J) 

'•IP.  I  TET  S  i  121  )-  TCOST  »  COST 

121  FORMAT (./Z.7K  TCOST *-••*  lPElg*.P  4 AX«pH  COST**  1PE1S.8//) 

2  o  ONCOST -TCOST  )  /DT.COST 
I  F(  0  »'0  o  GE  .Ho  OP..  K  *Gc*  >V.<)  H=XK 

A 

/  — 

FINDING  I*iTF.RVAL  WHERE  THE  DERI  VAT  1X/E  OF  THE  COST  PASSES 
THROUGH  ZERO 


DO  <*  :l-  i  4  MSTEP1 

4  COMTk  (  1  )  sCOVTRT I .)  -X:<---S'M  )  - . 

DO  6  T  k*-1  4.20 

DO  20  K---1WNSTEP-1' 

20  COS'TRT  1<  )  sOONTR  ( K )  ~2  »0**  (  I  - 1  >  vh’-^S  ( X  ) 

WRITE  (6-1  1*05  >.  I 

1-05  FORMAT (22H  TWO  IS  RAISED  TO  THE  -4-13-49H  -1  POWER) 
CALL  1 NTEGTT I ME 4  X, CCNTR ) 

CALL.  G'V ^NTtT  i.SF.i  X  4  CONT R  4 0>!DU ) 

DQOST=O.C 

DO  S  X  =  1  4  MS.TEP  1 

5  DCOST-DCOSTvSTX )  -"-CHOU  <;o 

WR !  TEC  6  1  !  07.)  COST  4  TCOST  ' 

107  FORMA TT/SK  CO  j>T  =4  ! PE 15 • 3  4 7H  TCOST=  4  1.PE1 5*8/) 

DO  21  l<=  1  4.VSTEFI 

21  CCNTR < K  )  aCC.NTR <  K )  r>2»* 0** <  !~T)*H-*STXX 

I F  (DCOST «L7 1 0 .0 )  GO  T0;  7 

I F  (I  .  EQ.4-20. }-  WR I  TSC-6-.r  1 00) 

ICO'  FORM  ATT  1'8H  TPCUELE  IN  SEARCH.) 

IF  (  I  .EG  *,20  )  STOP 
TCOST=COST 

6  DTCGST  =  OCQST 

cubic  interpolation  of  ><  . 


c 

c 

c 


A~2«C-::-*U~i  )-'-H 
3=6*0 

I  FT  I  *NE-«  1  )  ~=2*  C-:h:  <  I-2  )*H 
NCI  =6 

Z=  DCOST +0tC05T*3 . 0-"-  (  COST-TCOST  >  /  ( 3- A  5 
V/ =  7.  *  # 2  r  D CO  ST  *0  “COST 

v;=  sort  ( w  j 

XK=.»- ( »-A  >  vTSTCOS f-H-.'-Z ) / < DTCOST~DC03T-r2.0*W ) 

=WR.IJ1E,(^4-1TLL)^,XX.*JL<JL^  -  _ -J. 


101  FORMAT  (  AM  X:<=  4  iPEl  5 *3  4  3X 4  3H  A= «  VPE 1-5*  8.4  3X 4  3W  3=  4  IPETE  tS//) 
DC  10  I  =  1,4  NSTf?.l 

10  COXTPT  I  )  -CO:\TR  (  I  )  -!- XK*S  LI  )  - 


% 


■ 
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IF  (NCi  *EQ*  1  5  GO  TO  T3 

STonr£=eosi 

C  ALL.  I  N'TEG  <  T  I  ME  4  X  i OONTR  ) 

costxk=cost 

COST = STORE 

'  '  -CALL  GRAp.VT  { 1 1  ME « X « CCNTR  ♦  OHDU  5 

I  F  (  COST  *  LT- » COSTXK  vAK’D  «  TCO.ST » LT  t  COST.XK )  GO  TO  tfg 
TDSL-P 

,  -  00  22  I  -  1  i NSTTP i 
22  CONTR (  I  5  =CONTR  n  1  -XKtf-S <  !  5 
>  -  DO  25  1  =  2.45 

FAC*  I 

•0EL'kB+PAC#(.A^S>/5«<0  • 

-00  26  K~  1  4  NSTEP-1 
26  CiON.Tf?  ( ;< ):=CO\'TO  (  K  )-rDEL*5  < K ) 

CALL  I-NTSG ( T i MS 4  XVCOWTR ) 

CALL  GRADNT < T I  MS  4  X « CONTR 4  OHDU ) 

*  .  0C0ST  =  0*0 
-  bo  2a  x=i  VNsfEP.f 

28  DeOST=OCOST^DHDU(K)5'rS(iO 
DO  29  K={ iNSTEP 1 

29  -  CONTR  (  !< 5  r.  CONTR  <  K )  rDEL*S  <  K ) 

I F  ( DCOST  »LT  •  6  »  C  o  OR  »  I  . f£Q « 5  5  GO  TO  30 
TDEL*DcL.- 
0TC0ST=0C05T 
25  TCOST=COST 

30  A=Oc£L  . 

BsTOEL 
NCI  =  1 
GO-  TO  9 

12.  COST =COSTX’< 

•RETURN  ‘  - 

,13  CALL  1 NTEG  { T4  ME  i  X  rCCMTR ) 

-  CALL  GRAONT  ( T I  ME  4  X  i.CONTR  4  DHDU  )■ 

RETURN 

END 

SISFTC  GR 

SUBROUTd  K'E  GRADNT  ( T I  .ME  4  X  4  CONTR  4  DHDU ) 

COMMON/f  00/  DT I  ME  4  NSTEP  4  NSTATE  4  S  AREA  4  COSTi  FT  I M  i.FCON 
DIMENS  I  ON  X.( 8  44  00.)  4 CONTR  <  100  >  •  T I  ME  <  1  00  5-4  DHDU  <  1  00  5 
NSTSP1 sNSTE^+ 1 
'  DO  1  1  -  1  4  NSTEP  1 

.  ;  TlMs TIME <1,5  ,  -  - 

IF<TIM«C-T.T,I  )  R=0«0  :> 

I F ( T I M • LE • T I  5  R=RO 
X3=X(.34  I  )  ‘ 

GRAVT Y=  f  «  4 1:0020 19E.16/  ( 20 •  9ECH-X3  5  _ 

XMASSsFMASS  ( T.-TM )  J  - 
‘  THRUST =FTHP.US  C  T I M ) 

-•  RHO=FRHO(X3)  • 

■VS=FVS(X3)  ' 

XiM=X(  1 4  I)  /VS 

»  ;CLA=FCLA.(XM)  -  * 


TR .XVOST4RG4TI 


GRADIENT  EQUATION 


u  b  u 
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DHOU  <■  I  )  =~X  <  5  4  1  )  <  THMJST*S  1 N  <  CON'T^C  I  ) 
H-X(o4  i  ) *  < G*CLA-'-THRUST*COS.( C0N17R  (I  1}  ) 

a~R-;-‘co.\'TR  <  n 

DO  2  1=1 »NSTEP1 » 4  . 

WR 1 7F (64  100)  T I  MS (1)4  DHDU (  I ) 

FORMAT ( F 1 0. 2  4  5X  4 1 PE15 . 8 ) 

RETURN! 

END 

SUBROUT  I  NLr  SQUAT  (  TX  4  DXD7  4  T 1 M  4 U ): 
COMMOM/TOO/  97 1  ME  4  NS7ED  4  NSTATE-4  BARE  A 
D  I  MENS  I  ON  TX  (20)4  DXDT  (  20.)  4  Y  ( 1 0  )  / 

XI  =TXM  >. 

X2=T X(2)  - 

X3=TX(3) 

X4=TX(4> 

X5=TX(5) 

X6=TX ( 6 ) 

X7  =  TX  (  7  )  -  ‘  -  -  -  .  . 

X8=TX(8) 

GRAVTY= 1 . 4 1 00201 9E  1:6/ <  2C  »9E6*X3)**2 
T  H  RUB  T =F  T  H  R  U  S  ( T I M ) 

XMASS=r MASS ( T I M ) 

RH0=FR!-!O(X3) 

VS““VS(  X3 )' 

xm=xi/vs:  •  '  , 

CLA=I;'CLA(XM> 

coo coo  (>:?•: ) 

-6=  3 ARE  A*RHO*X  i  **2/2:.  O' 

DRAG=  ( CD0*CUA#U*#2-)  *C 
XL  I  F  7=  GL  A  4:-U‘"*0 


J-+GL-A'* 

/  ( X  (-l  4 


CON'TRX 
I  )  rrXMA 


T)*2  4  0 
SS<)  ■ 


*0-1  /XMA< 


4  COST-4  “7  I!M  4  FCONTR , XVOST  ?0  4  T  I 


’F I  NO  1  N'G  THE  DERIVATIVES  OF  RHO4  VS  4  CLA  4 


y ;  *• 


DO  2  :<=l44 

I.F  {!<  .  SO  *2*  AND  ..4 .  E4  .'LT  •  X3  *  AND  . X3 . LT • 6 .  E4  >  DVSDH  =  0  •  0 
I F  ( K  tlEO  •  2  »  AND  « .) .  S4  *  LT » X3  •  AND ♦  X3  . LT  . 6  «E4 )  GO  TO  2 
1 F  ( :<  4  FO » 2  •  AND  «  1  -i  6k’E?  .  LT .  X3  *  AND  ..X3 .  LT*  1 . 7S5).  -DVSDH=0  •  C 
I F(;<*E0«2« AND.  1 .6E5.LT .X3«AND.X3.LTi  l  »7E5)  GO  TO  2 

IF^K.EQc-3*  AND.X-V..LT.  1 .0)  DCLADM=0?i:0 
I F  ( K  «  SO  •  3  •  ANTJ  •  XM  •  LT  •  1  •  C )  GO  TO  2  :  . 

I F ( K  *  SO  •  4-»  AND •  X‘-1:.  L-7  •  0  »‘6' >  DCDODM=O.G 
IF( !< . EG . 4 . AND . X.M.LT .0.6)  GC  TO  2 
IF  <!<»LE» 2 )'  DELTA* =750. 

I F  ( K  •.GE’*  3 1  DELTAX=0»06 
.00  1  I  =  1.4  9 
FACT  =  I 

I F  <  i<  »LE  »  2 1  XT  =X3+ DELT  AX*  <  “  ACT-5'.  0 ) 

IF  (K . GE. 3 )  XT=X?:+  0 ELT AX#  ( FACT-5 •  0 5 

i  f  ( .  ea .  1 )  v  <  i  j  =frho  (  xt  > 

IF(K.e6.£)  Y ( I) =FVS (XT ) 

IF  (K. F0.3)  Y  Ci  )=FCLA(XT  ) 
t F ( K.  EO » A )  YCi 3 =- CDC (XT  J 

.DVDX::Xr423«^y^5,12:AE^^iE3-33.j^(.2.)^I;6..d^(,SA^34^(5i;k)vl:(;: 
bYDX=DYDX-M-3  .  *W.4.1  - 1 0 .  Wfn  >  +  H  *  *-Y  (5  )  “6  .  #Y  <  7  >  + Y  CS  )  )  /  ( 
DYDX=OYCX-r  ( Y  ( 3  )  -5 .  *  Y  ( 4  )  *3 •  *Y  ( 6 )  -V  (7 )  )  /•(  1 2.  ^DELTAXJ 
OYDXsDYDXv (  - Y (2  ) *6 .  * Y ( 3 > ~  1 3 .-”Y (4  )  + 10. *Y < 5  ) *3 YT 6 )  ) •/.( 
DYDX=CYOX-r  ( 3  •  «*Y  Cl  1  -46  .  **Y  ( 2 1 H-36*  -”Y  ( 3 }  -4B  .  «Y  ( 4.1*25  i  ( 5 


;l:2.*JiDE 

ll.-KDE 


L7AX.L 
LTAX ) 


1 2  .-^DE 
)  )/ 


El.  TAX) 


t 


yV».  »■  J»*J  'W* 


.*«  uU>V~  n/rt.  MtflU. 


3C/se/?o^5. 


.i  cir  **d§l~'a>o 

0Y,QX*DYdX/5*0  .:.  .  . 

JFtXt'ZQ*  I)  DRMODHaDYDX 
I F  ( K ,-  !-:n.t  z  )  DVSOH  ".DVD X 
IPCC.FO.S)  -  0CI-ADP--DYDX 
S  -  I  F(  ,< .  EC,r-V>  1 3  C.  f  >  0  D  M = D  Y  D  X 

V* 

C.  TME-:  PARTIAL  DER I-VAT I  VPS  OF  LIFT  AND  DRAG  V/*P»t«  XI  AND  X3 
C, 

•DDDVsCfl  (DCD0D^:/V3^UA^g-voCLAD--./VS,)  v:-:!-;0i:Xi*iA^rA"DR/.C/C. 

T  DLVC>yrD«.(.CLA*u/XJ#*2+U’'DCLADM/('yr.*Xl  )  > 

DL0H=G-T  f  <  X 1  r.RHO )  •-u::-DC.LADM*DVSD!-i/V 3J< -2  ) 

-  PDDK~0^.('-.Xi->:'E)C30l5i  -SDVSGK/VS*  Xg-U* *2* CCLADMW?  VSDW X.hA/S  **20  •:■  • 

-  :  LDRACivOnHCDHXRHO 

_  -  - 

C  the  .aiFFERENf  1-AL .  EQUATIONS  OF  f HE  STATE  AND  .ADJOINS  -STATE  V& 
C  -  •  >  =  - - ••  -  -  ----  -  -  .  -------  -  -  - 

DXDT  (  i->.aTHRUST«-e03  ( U )  /XMASS-DRAG/XMA3S~GP.AVTYY-S  2  N  (  '/.?. ) 

DXDT {,&•)  -XLI  FT/  <X.MASS*X1  )  +TBRU3T*S 2 N  ( U 3  /  ( X  4AS'S*X lJ-G3AVTY*CpS ( X2  j 
1/Xi  ' 

\  D  XDT‘<  3 )  -  X  i  -X-S  I  N^X3  -)  - 

DXDTK )  =X.r-:,-CO^  X?.) 

DXDT  ( 5.)  -XS>DDDV/Xf  .ASS-/6*  < DLVDV/XMA3S^"H".UST-"-S  i  N ( U )  /( *V,ASS*-X 1* *2  > 
1+GPAV7Y*c6s  ( X2-)  /XI  **2}-X7*S  2  N  C-X2)  ~X8*C0S ( X2) 

DXD-T<  6)  =X5«-GRAVTY*G0S(X2)-X5;*GRAVTY»SIN(X2  )  /XI  ~X7*XT*C03  (X2  )  + 
1X0*X1*SIN(X2)  • 

-  *  DXOT  (7J--FX5f1DDt^W*-^G^DLD.H;/.XjP/XMA'S^  P 

PXQT-(fS)  aO:*0\-v 
RETURN-  ’•*  '■'£■- 

^ND‘  '  ,  r-;  O  V'-  ' 

SIP.FTC  -"MAS  -0T'£:<  „  - 

:V  -  FUNCTION  FyASSXTIM )  I:  ’  - 

'GRAVTY=32>J:v;;  '  •  -  -  \  .  . 

/~  0  -  ----  J.  _4  .  .  ; 

y .  _  _  - , 

t;  8ASS  AS  -&  ^FUNCTION  OF  TIMS'  AND-'  GRAVITY 

^  '•  :  -  -  '■  \  -  •  -  -  -.- 
v  "  c  -  ,3  - 

_  I r  ( 0 •  0 •  LE  ».T 1  M  .•AND'.* T I M  .LT* 19.44)  FMASS=  t  1734  ,  *-S  1  «,25*TT M  ;)  /GRAVTY 

IF'{vl 9 » 44 *LE*  ? I M  .AND. TIM  *LE»28.4>)  FMASS“7-37 •  7/GRAVTY 
1-<1  ♦6034.7E~3*TTM  ».*S-3  •  2t>?7E~  1  *7PM  **7+29  »  1 634*f  I.M  *-*6-4  *'4721.39 
223*7-  I‘M  **5+4 .61  362£4*T  I.M  **4~? •  1  9 16383* T  I M  **3+t.  13573 1  k  7*T  I'M  **2 
3-7 i 97260 1  E7-*t  I M  )  /  ( £40 . *GRAVTY >-  “1. » 0 1 2 764653.6 /GR A VT Y 

I F.< T IK "<•  GT. * 28  *Y« )  FMASS?S12*  10246/GRAVT-Y 
'  F  M  A  S  3 = F 4  A  S  S + T  2  6  6  ♦  /  G  R  A  V  f  Y 

.  v  -  c  U»S'N-  ^  ;  's  - 

'  '  2ND  V  /,  - 

SJoF-C  FFHRU  '-  .  1 

’  » t  FUNCT  ION  FTrlRUS  <  TIM)  -  - 

-  -  V  °  *  ’  -  * 

■  C  THRUST  .  AS  A  FUNCTION  OF  TIME 

c  '  <-  <  •  "  '  • 

IF  (T  IM  .LT.  19.44  )  FTHRU5=T*.23E4 - 

.1 F  { TIT *6=;  ♦  1;9  «  44  i  AND*  T  I  F: « L.5-.  30 T35  >  FTRRUS*  1*  £62752-2* T I M’  **7 


••£.  295345+T  ii-l  *-::-6-:-'r..7 49 S 0 6£2* T  I'M  **3-7*36CVp9923”  H  l A 
*71 M  *-*-2-2 * 75743926 *TI M  **2+2*27356227*1 1 M  -7*97760 15 


4.F,(T  IK  *'GT  ♦  30 .35 )  F<THRUS=  C\* p 
REWR& 


n  n  <d  i.o  o  o  o  to  o  n  n  k  o  o  o 
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:i  I  r-.F'TC  FRK 

FUNCTION  FriHriXS) 

RHO  AS  A  FUNCTION  OF  X3 

1 1-  (X3.LT  *0 .0  )  FPHOi-O  . BOO  i  53- -  i :  “X.  3  '  ~  *  “D0927E~ i  3*-X3**2 

1  “6 . 9 60 B 5  7  F  ~ .-;*. -V X 3 •!■  2  «  376S9EU?*-3 

I F  ( X3  iLTn  .77232:::?  )  FR:  !C~  c 2370“-?  {  1  »  0- #  '  ;  '-i~“vX3  )  tt*4 *236 

I F  ( X3  tE0«  c  37322;: 15  >  FR HO-  .  7277-7 
TEMP-  1  o  4  J2  -I-  CX3  -«.  33o;?2E5  ■  / « 90-7-  7' 

1F(X3  cGTe  «3E73£,T7>  F RHO"  (  1  *  02  "  *7  '*  '*  7 ~  -  5 /,rX“ i'TcMP ) 

RETURN 
END 

1 BFTC  FVOS 

FUNCTION  FVS  <  X3 ) 


VS  AS  A  FUNCTION  OF  X3 

I F  (X3.LT  *0*  9  )'  FV S=S . E88SS4 E— 1 3*X3--‘ 

1-3 4 64 ??,B  1  .“-3-;'  X2-:- 1  .  i  i  6 87 92 3 
I  F  (0o  «L5  *X3  « AND  jX3  *LEcA«F4)  FV? 

)  -2  *  37S202E- 1  2*X3  *T 3+2 *014  3937 -2 •‘•72  “  :P- "  *  03  760-4"-:v*X3 

2+1  *  1 1691 ir 3 

IF  ( 4 « Gc  >'!  * LT  « X3  •  AND « X3  cLT.6rOE4 )  “  V7  ~ 7 *  N44 7T2 


T )  .77  ‘  <  •"7F“‘72'.X3-::-x2 
^  { '■••/  ,o  i;7-»c: 


-  2  .•  '■.  "3  °  i  E  -  2  2  72  *#5 

-1H  .1  ..04  -'r>"'6“  •,-<7*.Xn  &«g 


IFC6.k4.LE.X3  .ANi:.eX3  .LTd^C)  “VS 
1  **. . « 932 1  -  ?F~  1 7-:-x3  ••>.-4+6  2v2:  l  02F-1  ;</ 

2-2  o  1  4 773<N->3*'<~  ■>• }  *  1  1 632922 

IF  (>C3  *ea«l‘.“S)  FVS  “1..004SE3 

I  F  C 1  « E3.LT  .X3  « AND.X3  .LE«2»E-:vi  CVS  -  1' e-7C0396E- K 2»X3,  *H 

1-1  i23S934i*-12*X3  * *3+3.20960  /E~7*\3  v^2~3.  1  7^?9E~24X3  + 

22  *057099!:  3 

TF(i,.6EE>«LT»X3  .ANlVXS  *L7  «  1  *~/££ )  FVS  - 1  «,  1.057F3 
I r  ( X3.»  GT « 2 •  E5 )  FVS”  1 032* 


RETURN 

END 

I BFTC  FCL 

FUNCV I  ON  FCL A ( XM ) 

CL A  AS  A  FUNCTION  OF  MACK  NUMBER 


•Vor  m 


IF(0.  .LE.XV.  .AND.  XX  »LE.l  vO.)  FCLA-l  o02r-i 
IFM.O.LT.X'M  .  Aw  ) .  XN  .LT.1.65  FCLA -~C « 73371  SE-iSSXX 
1  +  1  »  736429E—  1  *X*^  #  4631  77E-.3  +X-"  -  2+3 * 7.11  7S9E~2-*XM  - 

2+1 .025505E-1  - 

IF  (1.6. LE .XN  « AND < XN  «••_£. 6  ♦  0 )•  F CL. ■ : - 2 * 4 O l-tog KE- 4-*XM  v*'4 
1  +4 . 54r 963E~2*XM  +42-2.637706E-2 :>X.'.  ^ 2 977S8/v--2-+X.Vi  + 

29.9930 1 7E-? 

FCL A-FCLA457 . 2953 
RETURN 

ENO  ,  - 

TSFTC  FCD  ■  -  ; 

FUN'CT  I  ON  FCDO  ( XV. ) 

COO  AC  A  FUNCTION  OF  MACK  •  NUMBER 

IF  (0.C.L2.XM  .AND.XsM  .iLF.C  »6>  FCDO  - 1  TS2  - 

I F  ( 0 . 6  .  LT .  XM  .  AND  »  XM  »LK  .  5  «  1 5 )  FCDO  \  1 . 33  79-2-KXN  K*4-  - 


*>?« 
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1+33*  V:C7-::-XM  ”>U-\  UC 11 r>7--~*XM  —4  *'$9Z40'2r- 

'  -j.F*,(-:i:*4SVL.T«*XM>  FCD0*-5, *»!?«•  ir-r  .j 

'  _  l--7;‘..r5081 .6E~2^XM  **3+0. .  ."7s-:’5'  ■  X*-  *  'T-7  i  ?  i 0A'6~=t-  T*TO 

;  2+6*9904876  _  / 

’I  P  ('4  »<?72  *  Lt  » XM  5  F  COO- OV 1 9 . 

V  -  -  'RETURN 

SfSFfC  .lK‘Ti“  - 

^us.'jcurr^ 

COMMON/TOO/  OT I  *'.£ « KSTfXV irT.-.T;. ,  •  -‘“A ,  C C *?  T » FT !  M  *  FCONTR  *  XVOST « RQl*  7 1 
•'0 1  PENSION  X  (B  <  TOOT  *  CCNTTI  K3i  n  ‘"U  00  )  i  f- !  TO  *  A )  i  TX  <1  0  >  »  DXOT  ( TO  ) 

:  .  M-T  *• 

NEG~0  r  l''_  : 

■</. ::  ~-stcre~u6 

•NSfE»lJ=NSTEP+l  ‘ 

DO  14  1N3=0  2  -  - 

.DO  TO.  JJr-T-rK'lfEP;  * 

*  -4-p-(^DVfe^.rl-*-Ak,e-*-&T0RE;#L7«X)>*0*7  *  "30  TO  F=T  ' 

G  .  -  #iMfiOT4^G  RUNGE~KuTTA  FO~v;;LA  7  - - 

c:  X',(N+  i  ),  =  X((N  )■  +  -<kr.O  +  +  2TX-P  +  i<~3>/6 

c  .  '  -  '  - 

•i:F>DiNO,..EO-i'l  )  J~JJ  *- 

.  I F (I  NO  <  EO  *  1  )  J 1  =  J  J  +1 
.  i:F;i?fNC«EG:*.l  )  H=D7  TMc  _ 

!F(iN0*E0*2)  J=NST.EPT'*:{  JJ-1  !  .  ,  ' 

'  \  F-CJ  ND  •  EG  *  2  )  J  l*«N STEP  5  “  J  J  . 

dKF-  H-WD  •  EG  «  2->  H=TfME(  J1  )-7;F..E{  J’T  ^0> 

v&  Fvl.K'OTNG  „  K~0* 

c  =  .  O  -  ,  ,  "  -  % 

'  ;  .  F.T;1  ME-T  i  ME  ( J;)  -  .  :  •<?>  . 

:  DO  i  K= 1  INSTATE  c  ’  •-  ^  .  I  - 

T /  T:X4kl  ~X>{K.I  J  )  *  .  : 

-.CSN-fRCsCONTPi  J)  T  5  -  ~ 

C  ALL  EGUA7  ( TX  fDXpT  |FT  TME  *  C0N7&0J  * 

do  2*  K=M.iNS:T*ATE  -  ;  -  '  _  '  *  -  :  . 

iZ  "  c  ( Kiji  )  -  DX  DT  <  :<  ),-”•  H  "  "  >  C.\ 

.P--  .  ■  ■  ■  .  ■  -  -  - 

C  ,,  FINDING  !<-l  ,v 

•!C 

:  ..  -F-T.I  WS«T«k'E  { J-)  +0  .-5*H  -  ;  ” 

I.F-t  I‘ND « F:G  *1; )'  J?=J+.l.  ,  *  . 

.6l  .iF(tNDic6.2>.  JZ=J-1 ^ 

tONIRps CCN7S'.<\'2;)  +-J FTIFF-T I M”{.J2 )  )  "  ( C.ONT-P  <  J ) -GON~P (.02  )  )  / ( T IKE  (  J.). 
L^TlMii.(;J218  - 

.  DO  3  k^l'iNSTATT 

3  TX  <  K1  =X'C<  *  J  )  +0  *  5  ::-F  ( :<  *  i  }  : 

i  ,  Cll^^l^AT^^Xfl^DT^FTI^E^GNiTJGLT:  ----- 


1  ,* 

4 

F(X.»  21  =0X07- CO 

I’- 

>C 

-  'C-,- : 

>1F  18  Oil  is'  GJ_JK  « 

i 

6 

«■> 

DO  f>  k=ilNS7A7 

■“  tx  <k )- X(i<4  41+6 

~  ' 

CALL  EGUAT ( TX  < 

«■  -  ; 

i 


fcc 


fu 


i  : 


•*  ,  1 


it- 


;*4 


B' 
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DC  6  |<=M  4  NST  ATE 
F  { K  4  3  )  *:0>,TT  < K )  #H 


C 

C 

C 


F  1  N'D  I  f  vo 


;c 

G 

c 


F  T I ME = T  IMF  (  J  )  +  H 
I  F  (  I  ND  *  HO  *1  )  CONTRO=CONTR  <  J6  1  ) 

.1 P  {  l'NP*  FO  »  2  )  cetJT.po=CON*TR  (' J-  1  ) 

DO  7  ,K-  It NST ATE 

tx  ( ic  >%=x  ( !< « J  >  -:-F !  k  ♦  3 ), 

CALL  EQjJAT  (  T-X.4 DXD~  ;*FT I  ME 4  CCN7RC  5 
DO  3  K=y  4  NST ATE  /  : 

F(.:<4  4-)=dxdt(<)  "*m 

FINDING  X ( Nvl > 


f  j 


.  DO  9  K?M INSTATE* 

9  X<l<.rJ4:).“X  ( K,»;J.).-5-.(.F .(.K.4.1. ).-h2»-p" F, (.!<,». 3 1+^'O-K”  (IK  4.31  *' ?F.VKs£M.j,6  ».0. 

S'TCDF=X-(3*D1) 

1C  .  CONTINUE 

GO  TO  (  1,5.4.14:)  4  1  NO 

25  DO  26  1= 1 4N3TATH 
CvDO  26  J=JJ  4,NSTFP 

26  X(I« J-XM-XCIm J)  ( 


C 

C 

e 


EVALUATING  the  approximate  final,  time 

15  DO  16  J,la2i;NSf.rPl  "  V) 

1  F  (  X  (  L’J  1  )  «LE  «  0  i  0 )  GO  f.C  17 

16  CONTINUE  .  -v  \T  “ 

WRlTFCf  4-100)  ’  ,  \  •  - 

ICO  FCRMATX//29H  /INCREASE  TrME  OF  I'NIECL  AT  I  ON//) 
TN3=NSTEP  ’  :  ‘ 

FTTM~-TNS{:-PT!.ME  . 

C0ST=X(*4'4N£TEP1  j  ‘ 

F  c  ONTD  =  CCNT  R  ( NS7EP  14- 
GO  TO  22 

17  FT  I  M=  7-l.Mc.-jt  J  l ’)  -  (  T.I  ME  (.  J 1  - 1 .)  -7  i'.ME  ( .  f \  5  )  *X(  3  4  J  l  )  / 

1  (X(34  Jl~4  )~X(3i  J1  )  )'  . 


C 

C 


AcPROX I MAT. NG 


COST  AND  CONTROL  AT  FINAL  TIME 


1=4 

CC ST  =X  <  1 1  J 1  -1  )  +  ( X  U  4  J  l't  1  )  -X  ( T  4  J i  “  1  )  ),'F( i- 1 1 M— Ti  MEX- J 1  - 1  >  )  / 

HTrME(Jl.U)  “  TIME  (J  1-1)) 

F'COWTR  =CCN7R ( J 1  r  1)  +  ( CONTRX  J 1 61  )  -CONTR  (  J 1  - 1)  ) •* ( FT  I M-T  I  ME<  J 1  -  iiVli/ 

1  ( 7 1  M3  (  J  1 6 1  „•)  6  7 !  MET  J  1**1  )  ) 

i  (ji-i  >  > 

FCONT  A- = CONTR  (;  J-l  - 1  )  6  C  GCN7R  ( J 1  )  -COA’TRX  J 1  - 1  )  />  *-(,FT  I M-T  1  M~  (  J 1  - 1  )  j  /  (IT TME 
CONTR  <  U1  )!=FCO\'TR  : 

NSTEP1=JI  :-r-  ' 

NSTEP  =  NSTE~l,r:l  -  - 
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12  'XtUNr'-^PlV-Xf  I  «N3J  - 

1/;  GCNTi  Nt,'-.  '•  '  -  .  "  . 

s '  NStAT{£-‘K'3  TAT  E/2 

'• .  XVOST^PR  I CE  ( CONTR-) 

Dp  PT= 0  *  5-::*G03 1  -”-*2  -XVOS  T 
‘-/-.DO  30  I  -  1  ♦MSTEPi  » 4 

30  \to  I  YE  (6<103)  T I  ME  (  I  >  4  X  ( i  ♦  I  )  (  X  {  ?  *  I  )  ( X  ( 3  *  J  >  <  X  ( 4  (  !  )  ,  CO\'TR  (  1  ) 

1  (  X  (5>«  1  )  ( X ( 3  (  I  )  i  X  ( V  (  /  )  ( X ( r  t  I  5 
1:03'  5F ORMAT ( 1  0  ( E  i'O  3 ♦  2X )  ) 

RE  TURK' 

END 

1‘SFT.C  PRi'C  DECK 

r  ONC.f.'l  ON  PRI.CF.  (CONTR) 

COMMOk'XTOO/  DT  IMS ,  NS7EP ♦ MPT  AT’  ♦  3  '.REA  ,0037 « FT  I M  ♦  FCONTR  (  XVOST » RO  *  T I 
DIMENSI  ON  GONTRMOG-) 

N=  7 1  /D  7 1'ME-'r  3  *  0 

RaRO  ,  '  t  *  : 

PR!  CEaO  o STRtfCCMT  R.<  1  )**2 
PR  TCEsPR  1 2  *0 "R-::-CT)T.!.7RT(  2:) 

DO:  1  1-3*  N_(  2 

FAC =1-1 

71  ,V,=~  AC-”-D“  I  ME 

I F ( 7 1 M • GT  «  T I  ).  R=0  »  0 

IF(T1M.LE»-TI,)  RsR‘0 

PR-1  CE=PR  1  CF+R^CON.TR  (  1  )  **£ 

TI  M=T.l  M-rDT  IsMS 
IF  (7IM*C-T.TI)  R=0  4  0 
IF(TIM.LE.Ti)  R=RO- 
1  PR !  CEsP R‘I  CE+2 . 04R/CON7R  (  !-.•:•  1  ->  **2 

PR  ICE=PR!CE'»0t:IME/3 .0 
RETURN 
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Appendix  B 

The  Drag  and  Lift  Coefficient  Ec 


,  \  p;.l 
<C\r 
rm 

\ 

+  ''i 

s 

wr/ 

In  Fig.  IB^Xjj  and  y^  are 
the  body  axes  of  the  missile. 

L  and  N  are  the  lift  and  nor- 


/  J  tively;  E$  and  are  parasite 

_ y.  and  induced  drag,  respectively. 

Fig.  IB,  The  Lift  and  All  other  drag  components  are 

Drag  Forces  Acting  on  the 

Missile  assumed  negligible. 

From  the  diagram  it  can  be  seen  that 

p‘=  F  , 

(IB) 

l  ~  F  C0S(et) 


(IB) 


If  oL  is  small 

P;  ^  F  oL 

i  -  F 

and  therefore, 

p4*  s:  L* 

The  total  drag  D  can  be  v;ritten  as 

V  =  Po  +  Pi 

Substituting  Eq  (?.B)  into  Eq  (3B) 

p  =  Pt  +  L<< 

The  hohdimensional  form  of  Eq  (4B)  is 

c0  =  Cpe  +  CL< 


(2B) 


(3B) 


(^B) 


The  Hi'ft  eoeff icient  can  be  written  as 

et~  Ct  * 

where  CT  ,  is  the  trimmed  lift  curve  slope  at  a  given 

L  t  C-i. 

Mach  number.  Substituting  Eq  (6b)  into  Sq  (5B), 

ct>  -  %  +  tb*-  (7B; 

Mote  that  the  use  of  an  aerodynamic  moment  equation  is  not 
necessary, since  only  trimmed  flight  conditions  are  used  to 
Oefine  the  lift  and  drag  coefficients. 
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Appendix  C 


L'agrangian  Dt  f  f crentlat ion  Formula 


Assume  the  function  F(x)  can  be  represented  by  the 
Lagrange  interpolation  polynomial  (Ref  6:  48) 

£  4  pJ  4  +  BPl  < 

,k~0 


where 


■  ,>4  /  v  :  / 

lM  -  ■(*-&>'{*-*■) - -fr-aVi) 


f(*)  =  f *-£>(*- *4  -  •  •  •  {<*-  %) 


4  -  f -i  *k) 

yjfi  =  ii u/ntsr  fois/ts  ufej  iy  fAe  ApffOx 


f  (xo§  x^»  •  • .  .r j. ,  x)  is  the  divided  difference  of  the 

<k+1)th  <***•  f(^-  /fe) 

=  --y^r" 

Hwi  =  MSzMmI  . 

V  *  -  *» 


th*r~**M  =  - 

The  Xi*s  (i=0,  l,***k).  dre  spaced  equally  apart. 


*3  -  44  v "  **0 
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By  differentiating  Eq  (1C)  r  times*  the  r  derivative 
of  f(x)  can  be  approximated  as  .  . 

J -jfte  4  *  f'fe  „ 

/  A  _  K  T  0 

Assuming  E'  ' ( x )  is  small  and  therefore  can  be  neglected, 

Eq  t%C)  becomes 

t<%  --  i  ft 

tfhich  is  the  general  form  of  the  Lagrangian  differentiation 
formula  with -the  error  term,  neglected.  To  obtain  the.  five 
formulas  used  in  this  paper,  set  r  =~i,  n  =  4,  and  x  =  x^ 

(i  =  0,1,... 4)  in  Eq  (3C).  These  formulas  are  located  in 
the  wE4uatM  computer  subprogram  (Appendix  A) . 


Appendix  D 


Polynomial  Least  Squares  Curve  Pitting 


The  Legendre's  principle  of  least  squares:  given 
exact  or  equally  reliable  data*  assume  that  the  best 

approximation  of  curve  fitting  is  one  for  which  the  aggre¬ 
gate  of  the  squared  error  over  the  entire  domain  is  least 
(Ref  5:  63). 

The  exact  values  of  f(x)  are  known  at  discrete  points 
corresponding  to  xq,  x^,  . . . .xm  over  the  interval  (xQ»  xm). 
It  is  desired  to  approximate  f (x)  in  the  form 

•?(*;)=  (1=6,1,-*); 

k*0 


Define  the  error  r(x)  as 


Hh)  -  f(*i)  "  i 


k*0 


If  the  a's  are  determined  such  that 


it*  <3D> 

4-  Q  l~0  teO- 


X'S  minimum,  then  the  best  approximation  in  the  least 
squares  sense  Is  obtained.  The  minimum  of  R  can  be 
obtained  by  ordinary  calculus: 


GGC/EE/70-5 


ff-  Tr-fl  ^d=  °;  &hv~ rt) 

^  ^  La~-o  (5d) 

Eq  C5B)  represents  the  normal  equations.  Substituting  the 
a*s  obtained  from  Eq  (5D)  into  Eq  (10)  yields  the  fitted 
polynomial. 

To  increase  the  accuracy  of  the  curve  fit,  the  range, 
of  the  independent ' variable  is  normalized  over  -1  ^  1 

by  applying "Eq  (6D) , 

/  zxi-  y,n  ~  r-° 


- 


The  following  equation  yields  the  coefficients  associated 
with  the  unnormalized  independent  variable  (Ref  51  69)-. 

h  =  4*  (T)  t4  k*f  A  <7! 


where 


n  ~ 


c/eyee  of  ike  po\ymrAir/fi\ 

I 

Xrf:Xo 

Yi  jn+d  - :  ~  ‘(a-a  +  i) 

4  *  * *  *  *  2  *  i 


(?)-.*;  0 


r* ;  (!)  * ' ;  (!)  ~ 1 


==  +* 
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